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1 .  SUMMARY  OF  OBJECTIVES  &  PROGRESS 

The  objective  of  this  project  is  to  use  a  combined  Direct  Numerical  Simulation  (DNS) 
and  Large-Eddy  Simulation  (LES)  campaign  to  address  the  sub-grid-scale  (SGS) 
parameterization  of  stratified  turbulence  and  surface  prescriptions.  The  work  has  three 
distinct  components:  conducting  and  analyzing  the  DNS  solutions;  conducting  and 
analyzing  the  LES  solutions;  and  improving  the  SGS  description  used  in  the  LES  so  that 
LES  predictive  ability  is  enhanced. 

During  the  last  12  months  we  have  conducted  focused  analysis  of  our  DNS  solutions, 
explored  new  gravity-wave  breaking  simulations  (including  quantifying  the  turbulent 
Prandtl  number),  evaluated  the  numerics  used  in  our  LES  code,  conducted  new  work  on 
optimal  perturbation  of  stratified  shear  flows,  extended  verification  comparisons  between 
our  DNS  solutions  and  recent  aircraft  measurements  (and  found  excellent  agreement), 
and  developed  a  Bayesian  hierarchical  methodology  for  deducing  unresolved  turbulent 
motions  in  mesoscale-model  simulations  of  stratified  flows.  In  addition,  we  also  helped 
write  a  successful  DoD  High  Performance  Computing  Modernization  Program  (HPCMP) 
proposal  for  challenge  status  for  the  computer  work  associated  with  this  project  and  a 
successful  DoD  Golden  Opportunity  Capability  Applications  Projects  (CAP)  proposal  that 
will  provide  an  additional  1  to  2  million  hours  of  computer  time  on  newly  delivered  MPP 
systems  before  they  are  made  generally  available  to  DoD  users.  The  CAP  program 
benefits  the  DoD  HPCMP  by  having  new  systems  rigorously  tested  by  experienced 
users. 

This  report  focuses  primarily  on  significant  increase  in  performance  of  the  new  LES  code 
we  have  developed  for  this  project.  We  gauge  the  LES  performance  through 
comparison  with  our  DNS.  This  has  been  the  focus  of  our  efforts  during  the  last  three 
months.  For  more  detailed  discussion  of  other  achievements  occurring  earlier  in  the 
fiscal  year,  please  refer  to  the  quarterly  reports  covering  that  previous  time  period. 

2.  RECENT  ACHIEVEMENTS 

The  most  noteworthy  achievement  during  the  last  three  months  is  the  significantly 
improved  accuracy  we  have  obtained  with  the  LES  code  developed  under  this  project  for 
stably  stratified  shear  turbulence.  The  code  employs  Germano’s  identity  (Germane,  et 
al.,  1991 ,  Phys.  Fluids,  3, 1760-1 765)  with  two  filter  widths  (A  and  2A)  to  implement  a 
dynamic  Smagorinsky  SGS  model  for  momentum  flux  and  a  similarly  defined  model  for 
turbulent  heat  flux.  Specifically,  our  model  employs  the  following  definition  for  the  SGS 
vector  heat  flux 

where  C,  is  the  model  coefficient,  A  is  the  filter  width  for  resolved-scale  features,  \S\  is 
the  resolved-scale  strain  rate,  and  ^T/^x^  is  the  resolved-scale  temperature  gradient. 
Equation  (1)  is  also  evaluated  for  results  filtered  by  2A;  this  changes  the  values  for  q^, 

|5|,  and  dT/dx^ ,  but  It  is  assumed  that  C,  is  independent  of  the  filter  width  (this  is  the 

nature  of  Germano’s  identity).  It  is  worth  noting  that  the  form  used  in  Equation  (1)  is 
somewhat  arbitrary.  It  is  suggested  by  analogy  with  the  Smagorinsky  momentum-flux 
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Figure  1.  Turbulent  Prandti  number  versus  time  for  TKE  (dotted)  and  dynamic  (solid)  SGS  procedures  for 
Re=2000.  Filtered  DNS  data  are  shown  with  solid  symbols.  The  LES  runs  were  conducted  with  a  factor  of 
6  fewer  grid  points  in  each  special  direction  and  roughly  1300  times  less  computer  time  than  the  DNS, 
which  was  run  with  720  x  240  x  1440  spectral  modes. 


model,  but  it  is  by  no  means  the  only  form  one  might  select.  Therefore,  It  should  be 
understood  that  we  employ  Equation  (1 )  as  a  first  attempt,  and  we  suggest  a  more 
complete  analysis  should  be  done  to  examine  discrepancies  that  may  result  when  using 
Equation  (1).  Nevertheless,  as  we  demonstrate  with  Figures  1-3  below,  results  obtained 
with  Equation  (1 )  represent  a  significant  improvement  over  the  state-of-the-art  TKE 
model  with  which  we  began  this  project. 

The  advantage  to  modeling  the  flux  for  both  SGS  momentum  and  heat,  as  we  do  here, 
is  the  resultant  ability  to  compute  the  turbulent  Prandti  number  (Prturb=ve/K©,  where  Ve 
and  Ke  are  the  eddy  viscosity  and  eddy  diffusivity)  from  the  solution,  rather  than  being 
forced  to  assume  a  value  for  It,  as  has  typically  been  the  case  in  work  reported  in  the 
literature  (and  is  also  the  case  for  our  TKE  model).  We  have  employed  both  methods 
(and  reported  the  results  in  our  26  Jan  2004  through  26  April  2004  quarterly  report); 
however,  our  previous  results  suffered  from  an  implementation  bug,  which  we  have 
since  fixed.  The  new  results  are  presented  in  Figure  1 .  The  data  points  represent 
filtered  DNS  solutions.  The  solid  (dotted)  curves  result  from  the  dynamic  (TKE)  model. 
Prturb  for  the  TKE  model  is  a  specified  fuction  of  the  local  Richardson  number.  Its 
relatively  good  agreement  with  the  DNS  indicates  the  quality  of  the  calibration  used 
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Transport  coefficienls,  Dynamic  model,  no  upwinding,  120X40X240 


Figure  2.  Eddy  coefficients  for  momentum  and  heat  for  a)  TKE  and  b)  dynamic  S6S  models,  (k^)  is 
shown  in  red  (blue).  Data  points  indicate  DNS  results  and  are  identical  in  both  panels.  Note  the  dramatic 
difference  in  the  vertical  scale  required  as  the  TKE  results  grossly  over  estimate  the  eddy  coefficients. 


to  specify  Prturb  9  priori.  The  results  from  the  dynamic  model  are  in  excellent  agreement, 
especially  when  one  considers  that  they  are  in  no  way  constrained,  but  instead  result 
solely  from  Equation  (1),  the  Smagorinsky  model  for  SGS  momentum,  and  the  dynamic 
procedure. 

Figure  1 ,  by  itself,  would  suggest  the  two  procedures  (TKE  and  dynamic  modeling)  are 
roughly  equally  accurate.  However,  when  we  examine  v©  and  Ke  separately  (as  opposed 
to  their  ratio,  Prtuit)  we  see  that  the  dynamic  model  is  significantly  more  accurate  than 
the  TKE  model.  This  is  demonstrated  in  Figure  2,  which  shows  the  time  histories  of  the 
eddy  coefficients  for  the  two  procedures.  Figure  2a  shows  the  results  from  the  TKE 
procedure.  The  eddy  coefficients  are  grossly  over  predicted  (by  at  least  a  factor  of  5), 
indicating  numerical  solutions  which  are  far  too  viscous  and  diffusive  for  the  flow 
simulated.  Figure  2b,  which  shows  the  results  for  the  dynamic  procedure,  depicts  a 
much  more  accurate  solution,  with  the  SGS  quantities  slightly  under  predicted  (more  is 
said  about  this  result  below  —  it  is  not  a  general  result,  and  it  indicates  room  for 
improvement). 

In  addition  to  the  eddy  coefficients,  we  have  also  examined  the  2"‘^-order  structure- 
function  fits  predicted  by  the  LES  models  and  compared  them  with  DNS  results.  A  plot 
including  DNS  results  and  those  from  the  dynamic  procedure  is  shown  in  Figure  3;  it 
depicts  profiles  with  height  at  7  distinct  times  for  the  exponent  a,  coefficient  Cj^,  and 
inner  scale  £o  defined  by  fits  to  the  2"‘^-order  structure  function  for  temperature  of  the 
form  AfT  =  Cj  r“.  The  DNS  results  were  previously  validated  against  atmospheric 
measurements  (see  Werne  &  Fritts,  2000,  DoD  HPC  User  Group  Conference 
proceedings,  Albuquerque,  NM  and  our  1  Feb  2004  through  30  April  2004  report.) 
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Figure  3.  Structure-function-fit  parameters  versus  height  for  LES  (biue)  and  DNS  (red)  solutions. 
Structure  functions  are  computed  separately  for  streamwise  (x)  and  spanwise  (y)  spatial  separations. 
The  two  solutions  are  statistically  indistinguishable,  demonstrating  that  the  LES  is  able  to  reproduce  the 
DNS  results  while  using  1300  times  less  computer  time.  See  Weme  &  Fritts  2000  for  details  on  the 
computation  of  the  -order  structure-function-fit  parameters. 


The  results  of  Figures  2  and  3  indicate  that  a  factor-of-6  reduction  in  the  number  of 
spectral  modes  used  in  each  spatial  direction  with  the  dynamic  LES  procedure  can 
reproduce  the  DNS  results.  Since  the  time  step  scales  inversely  with  the  number  of 
modes  in  a  spatial  direction,  this  is  equivalent  to  a  factor  of  6'^«1300  savings  in  computer 
time.  Hence,  this  allows  us  to  now  conduct  LES  simulations  that  both  1 )  span  greater 
regions  of  parameter  space  and  2)  attack  problems  involving  atmospheric  dynamics  and 
physics  occurring  on  larger  scales  of  time  and  motion  than  is  possible  with  the  DNS. 

The  results  also  indicate  that  we  can  indeed  realize  our  initial  goal  of  bootstrapping  our 
way  to  improved  parameterization  schemes  by  employing  DNS,  LES  and  meso-scale 
models,  as  we  have  begun  here  (see  the  program  plan  we  laid  out  in  our  original 
proposal  for  this  work). 

When  we  initiated  this  work,  we  did  not  know  exactly  how  we  would  extend  results  from 
DNS  and  LES  scales  all  the  way  to  the  much  larger  scales  associated  with  mesoscale 
models.  Recently,  however,  we  developed  a  Bayesian  Hierarchical  Modeling  (BHM) 
approach  to  address  this  need  (see  our  1  Nov  2003  through  31  Jan  2004  quarterly 
report).  The  model  constructs  the  statistical  likelihood  of  mesoscale  SGS  quantities 
based  on  combined  computed  and  empirically  obtained  probability  distributions.  Such  a 
probabilistic  approach  is  the  only  manner  in  which  an  SGS  model  may  be  built  when  the 
resolving  ability  of  the  (mesoscale)  simulation  is  near  or  excludes  the  peak  in  the  energy 
spectrum,  which  is  the  case  for  stratified  turbulent  motions  in  any  mesoscale  model 
(even  high-resolution  models;  see  our  1  Nov  2003  through  31  Jan  2004  report).  The 
difficulty  with  such  a  BHM  approach,  however,  is  establishing  meaningful  confidence  in 
the  procedures  used  to  construct  the  likelihood  and  so-called  a  priori  probability 
distributions.  Nevertheless,  because  we  have  validated  our  DNS  against  atmospheric 


4 


Figure  4.  Eddy  viscosity  arid  eddy  diffusivity  for  LES  solutions  (open  symbols)  versus  numerical  resolution. 
The  original  DNS  results  were  computed  with  720  x  240  x  1440  spectral  modes.  D  scales  as  1/nx.  LES 
results  are  compared  with  truncated  DNS  results  (solid  symbols).  The  lower  panel  shows  compensated 
plots  which  are  multiplied  by  nx^  (to  remove  the  7^  dependence  of  the  eddy  coefficients  -  see  Equation  1). 
The  remainder  in  the  lower  panel  depends  on  the  SGS  strain  rate  S|.  Note  the  different  systematic 
dependences  of  the  eddy  coefficients  for  LES  and  DNS  solutions.  These  differences  between  DNS  and  LES 
are  understandable  and  acceptable  at  large  values  of  nx,  as  the  SGS  procedure  is  not  designed  to  work  in 
the  dissipation  range.  However,  the  diverging  results  at  low  nx  is  of  concern  and  warrants  further  study. 

See  text  for  additional  discussion  of  this  point  and  our  plans  for  rectifying  this  situation. 


observations  and  then  our  LES  against  our  DNS  (see  Figures  2  and  3),  we  can  have 
absolute  confidence  in  the  probability  distributions  we  constaict  for  the  subgrid-scale 
[uiUj)  and  (M,e)  momentum  and  heat  fluxes  for  stratified  wind-shear  events.  We  must 

now  continue  our  program  for  wave-breaking  and  for  coupling  processes  (likely  mitigated 
by  waves)  between  separate  atmospheric  mixing  layers.  However,  with  further 
refinement  to  our  LES  code,  such  work  can  be  more  complete,  can  be  carried  out  on 
much  larger  domains  (allowing  for  much  better  quantification  of  the  needed  a  priori 
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probability  distributions),  and  will  proceed  at  a  faster  pace  than  the  initial  DNS  studies 
conducted  so  far. 

Nonetheless,  in  order  to  realize  the  goal  of  accurate  and  reliable  SGS  procedures  for 
stratified  turbulence  at  mesoscales  which  account  for  ail  relevant  SGS  processes,  we 
must  have  confidence  that  the  statistics  derived  from  cataloged  LES  solutions  are 
realistic  for  each  of  the  phenomenon  studied.  Currently  we  have  achieved  this  for  wind 
shear.  In  order  to  obtain  similarly  validated  results  for  wave  breaking  and  for  coupling 
between  mixing  layers,  we  must  either;  1)  obtain  DNS  solutions  for  these  other  problems 
or  we  must  2)  demonstrate  that  the  LES  solution  procedure  we  have  implemented  does 
not  depend  on  the  phenomenon  studied  or  the  filter  scale  A.  Unfortunately,  our 
solutions  exhibit  slight  but  noteworthy  systematic  dependence  on  A;  see  Figure  4. 

Hence,  before  we  can  make  further  progress  on  a  general  SGS  formulation  for  stratified 
dynamics,  we  must  first  Improve  our  dynamic  SGS  scheme  so  that  the  A-dependence  is 
removed.  Currently  we  anticipate  two  potentially  promising  approaches  for 
improvement.  One  involves  examination  of  additional  terms  to  augment  Equation  (1). 
The  other  involves  extending  the  dynamic  procedure  to  include  a  third  filter  width  (at  3A 
perhaps)  so  that  we  may  evaluate  differences  In  the  spectral  slope  at  the  different  filter 
scales  so  that  we  are  not  forced  to  assume  scale  similarity  at  the  cut-off  and  test-filter 
widths  (as  we  do  now).  We  will  explore  these  alternate  approaches  if  funding  permits. 

We  have  continued  with  the  wave-breaking  solutions  begun  earlier.  We  now  have  a 
sufficient  number  of  runs  to  serve  as  a  database  against  which  new  LES  results  can  be 
compared. 

3.  PERSONNEL 

The  graduate  student  Matt  Tearle,  who  was  working  on  optimal  perturbation  problems 

related  to  this  work,  has  obtained  his  Ph.D.  and  is  no  longer  working  on  this  project.  To 
conserve  resources.  Dr.  Joe  Weme  charged  no  time.  Dr.  Tom  Lund  charged  only  56 
hours  total  to  this  project  during  this  reporting  period. 

4.  TRAVEL 

During  the  last  three  months.  Dr.  Tom  Lund  traveled  to  Williamsburg,  VA  In  June  to 
attend  the  annual  DoD  HPCMP  User  Group  Meeting. 

5.  DIFFICULTIES 

During  the  most  recent  three  months  of  this  reporting  period  we  have  throttled  back  our 
efforts  significantly  in  order  to  conserve  the  limited  funds  that  remain  for  use  later  in  the 
year.  At  this  point  in  the  project,  available  funding  is  significantly  below  that  anticipated 
by  the  original  contract,  so  we  are  forced  to  reduce  the  work  we  do  and  concentrate  on 
those  Items  we  feel  are  the  most  important.  This  Is  becoming  especially  foistrating  now 
that  our  dynamic  SGS  model  is  performing  well  enough  to  allow  us  to  perform  much 
more  realistic  and  Informative  parameter  studies  than  we  have  been  able  to  before.  This 
Is  exacerbated  by  the  knowledge  that  we  have  succeeded  in  securing  possibly  as  much 
as  two  million  more  hours  of  supercomputer  time  through  the  HPCMO’s  new  Golden 
Opportunity  Capability  Applications  Projects  (CAP)  program,  to  which  we  applied  during 
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this  reporting  period.  The  CAP  computer  time  is  in  addition  to  the  HPCMO  Challenge 
Award  we  won  with  others.  It  represents  a  significant  increase  in  resources  for  this 
project,  but  it  also  requires  significant  effort  to  port  the  code  to  the  new  IBM  platform, 
conduct  scalability  tests,  and  complete  a  series  of  very  large  (i.e.,  1000+  PE’s) 
supercomputer  runs. 

6.  PLANS 

We  plan  to  complete  three  journal  articles  currently  in  preparation,  participate  in  the 
HPCMO  CAP  program  (porting,  testing,  and  running  our  Boussinesq  code  on  the  new 
IBM  system  (kraken)  at  NAVO),  and  complete  analysis  on  new  higher-Re  gravity-wave 
simulations  we  have  conducted.  We  hope  we  will  be  able  to  expand  our  parameter- 
space  studies  with  the  new  stable-layer  LES  code  we  have  developed,  but  that  seems 
unlikely  without  additional  funding.  The  new  LES  code  represents  the  most  promising 
development  in  this  work  to  date.  Combined  with  the  new  Bayesian  approach  we 
developed  as  part  of  this  project  for  meso-scale  model  parameterization,  results  from 
new  LES  series  could  have  a  dramatic  impact  on  the  accuracy  and  reliability  of  SGS 
methods  for  Atmospheric  Decision  Aid  modeling. 
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